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Abstract. - A new approach to exploring low-temperature excitations in finite-dimensional 
lattice spin glasses is proposed. By focusing on bond-diluted lattices just above the percolation 
threshold, large system sizes L can be obtained which lead to enhanced scaling regimes and 
more accurate exponents. Furthermore, this method in principle remains practical for any 
dimension, yielding exponents that so far have been elusive. This approach is demonstrated by 
determining the stiffness exponent for dimensions d = 3, d — 6 (the upper critical dimension), 
and d = 7. Key is the application of an exact reduction algorithm, which eliminates a large 
fraction of spins, so that the reduced lattices never exceed ~ 10 3 variables for sizes as large 
as L = 30 in d — 3, L = 9 in d = 6, or L = 8 in d = 7. Finite size scaling analysis gives 
2/3 = 0.24(1) for d = 3, significantly improving on previous work. The results for d = 6 and 
d — 7, ye — 1.1(1) and yt — 1.24(5), are entirely new and are compared with mean-field 
predictions made for d > 6. 



Introduction. - In this Letter we propose to study low-energy excitations in Edwards- 
Anderson spin glasses using bond-diluted lattices. All previous studies of such excitations have 
been hampered by less-than-adequate scaling behavior that can be obtained with the limited 
systems sizes L accessible on undiluted lattices. Consequently, important questions regarding 
the spin-glass state in finite dimensional systems remain unresolved [1,2]. 

A diluted lattice exhibits less frustration on short distances, and it would appear that 
the trade-off between potentially larger, yet less frustrated lattices should not improve the 
attainable scaling regime. Here, we find for the case of the ground-state stiffness [3-5] of 
± J spin glasses on diluted hyper-cubic lattices that instead scaling corrections diminish with 
the dilution. Hence, in combination with finite-size scaling, a dramatically increased scaling 
window is obtained. This improves the value for the stiffness exponent from the previously 
accepted value of at best 2/3 = 0.2(1) [6-8] to our new value of 2/3 = 0.24(1), which may 
make the difference between merely knowing that a positive 2/3 exists and potentially using 
exponents to distinguish theories [2]. 

As a further benefit, the relevant bond-density regime is located just above the percolation 
transition, at which the mean connectivity for diluted lattices barely exceeds unity, invariably, 
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for any dimension d. Accordingly, using exact graph reduction methods, we can study exci- 
tations in all dimensions with system sizes L that are beyond the scope of undiluted lattices. 
For instance, reducing lattices with up to 9 6 spins, we obtain an entirely new result for the 
upper critical dimension d — 6, ye = 1.1(1). Similarly, using lattices with up to 7 8 spins, 
we obtain yr = 1.24(5), which suggests that the stiffness exponent increases above the upper 
critical dimension. For results in d — 4 and 5, see Ref. [9]. 

On one hand, an increase of yd with d would appear to be consistent with replica theory 
(RSB) predictions [10]. At least for d > 6 bulk spins should dominate energy fluctuations 
induced at the boundary (see Eq. (20) in Ref. [10]), implying y d — dp, where \i is the exponent 
governing energy fluctuations in the (infinite-dimensional) Sherrington-Kirkpatrick model. On 
the other hand, Ref. [10] also conjectures fi = 1/4, i. e. y^ SB = 1.5 and yf SB = 1.75, well 
above the results presented here. 

Measuring Stiffness. Stiffness refers to the ability of a spin system in its ordered 
state to resist the nucleation of domains of overturned spins. The creation of such domains 
may entail an energetic penalty for forming an interface. The more ordered the state, the 
higher the energetic barrier and the more "stiff" the resistance. In turn, if the creation of 
such an interface is accomplished with no or reduced penalty for increasing linear size L of 
the domain, one may conclude that the system is in a high-temperature state and poses no 
resistance against fluctuations spreading through the system. Thus, the stiffness exponent y 
(often labeled 9) is a fundamental quantity to characterize the low-temperature state of a spin 
glass [3, 11, 12], and its value is frequently relied upon [2, 13]. For instance, Ref. [13] arrives at 
conclusions excluded by the more precise value presented here. 

The stiffness exponent is typically measured via the "defect" energy AE of an interface 
induced in a system of size L by swapping between periodic and anti-periodic boundary 
conditions along one spatial direction. In a spin glass the interface of such a growing domain 
can take advantage of already-frustrated bonds to grow at a reduced or even vanishing cost. 
The width a of the distribution P(AE) provides a measure for the energetic cost of growing 
a domain of overturned spins. To wit, 

a(AE) ~ L y , (1) 

where L here refers to the size of a system with an inverted boundary condition. Clearly, it 
must be y < d — 1, and a bound of y < (d — l)/2 has been proposed for spin glasses [11]. 
Particularly, ground states of systems with y < are unstable with respect to spontaneous 
fluctuations, which can grow at no cost, like in the case of the one-dimensional ferromagnet 
where y = d — 1 = 0. Such a system does not manage to attain an ordered state for any finite 
temperature. Conversely, a positive y at T = indicates a finite-temperature transition into 
an ordered regime while its value is a measure of the stability of the ordered state. 

Due to its importance, there have been many attempts to approximate stiffness exponents 
in finite-dimensional spin glasses [4-8,14-18], using transfer matrix, optimization, or renor- 
malization group techniques. It has been argued long ago that y < for d < 2 and y > for 
d > 3 [4,6]. Only recently, the stiffness exponent for d = 2, below the lower critical dimension, 
has been determined to considerable accuracy, yi = —0.282(2) [16, 17]. There has been little 
progress for d > 3, despite significant increases in computational power. In d = 3, the accepted 
value so far has been y 3 w 0.19 [6,7], although there have been investigations recently pointing 
to a larger value, such as 0.23 [8] or 0.27 [16]. All of these studies are based on fitting power- 
laws over exceedingly narrow scaling windows at relatively small system sizes, 4 < L < 10, 
causing large uncertainties. Our value in d — 3 is at the upper end of previous estimates and 
amazingly close to (but distinct from) the Migdal-Kadanoff prediction, j/g IK = 0.25546(3) [18]. 
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To understand the shortcomings of previous investigations, it is important to appreciate 
the complexity of the task: Most numerical studies are based on sampling the variance <j{AE) 
of the distribution of defect energies P(AE) obtained via inverted boundary conditions (or 
variants thereof [16]). Thus, for an Ising spin glass with periodic boundaries, an instance of 
fixed, random bonds is generated, its ground-state energy determined, then all bonds within 
a hyperplane have their sign reversed and the ground-state energy is determined again. The 
defect energy AE is the often-minute difference between those two energies. Many instances 
of a given size L have to be generated to sample the distribution of AE and its width a(AE) 
accurately. Finally, a(AE) has to be fitted to Eq. Q for a large and asymptotic range of L. 

Even small errors in the energy for either boundary condition, by way of their subtraction, 
may soon lead to extreme inaccuracies in P(AE). While for d < 2 efficient algorithms exist 
to determine ground state energies exactly, and large system sizes can be obtained [16,17], 
for d > 3 no such algorithm exists: The minimization problem is NP-hard [19] with the cost 
of exact algorithms rising faster than any power of L. Hence, the values quoted previously 
for y3 were either based on small systems, L < 4 [6], or on elaborate heuristic methods with 
L < 10 that lead to statistical and systematic errors [7,8]. 

Reduction Algorithm. - To overcome those limitations, we propose to increase system size 
L without increasing the optimization problem by considering reduced, bond-diluted lattices. 
As long as the bond density p is sufficiently above the bond-percolation threshold p c , the 
dominant cluster is effectively compact so that the asymptotic scaling behavior expressed 
in Eq. Q ~ and the stiffness exponent y - is independent of the bond density p [18,20]. 
Furthermore, we have developed a new, exact algorithm, that is capable of drastically reducing 
the size of sparsely connected spin glass systems, leaving a much reduced graph whose ground 
state can be approximated with great accuracy. 

We will describe the reduction algorithm in more detail elsewhere [21], including its ability 
to compute the entropy density and overlap for sparse spin glass systems (see also [18]). We 
focus here exclusively on the reduction rules for the energy at T = 0. These rules apply to 
general Ising spin glass Hamiltonians 



with any bond distribution P(J), discrete or continuous, on arbitrary sparse graphs. For 
convenience, we use a ± J distribution on bond-diluted hyper-cubic lattices here. 

The reductions effect both spins and bonds, eliminating recursively all zero-, one-, two-, 
and three-connected spins. These operations eliminate and add terms to the expression for 
the Hamiltonian in Eq. (j2J), but leave it form-invariant. Offsets in the energy along the way 
are accounted for by a variable H a , which is exact for a ground state configuration. 

Rule I: An isolated spin can be ignored entirely. 

Rule II: A one-connected spin i can be eliminated, since its state can always be chosen 
in accordance with its neighboring spin j to satisfy the bond J+j. For its energetically most 
favorable state we adjust H := H — \Jij \ and eliminate the term — Jjj Xi Xj from H. 

Rule III: A double bond, J>y and J^j, between two vertices i and j can be combined 

to a single bond by setting Jjj = + J^) or be eliminated entirely, if the resulting bond 
vanishes. This operation is very useful to lower the connectivity of i and j at least by one. 

Rule IV: Replacing a two-connected spin i between some spins 1 and 2, the graph obtains 
a new bond Ji^, and acquires an offset H Q := H D — AH, by rewriting in Eq. J5J 




(2) 



<i,3> 



Xi{Ji,\Xi + J lt2 x 2 ) < \ Ji,iX\ + Ji,iXi\ = J\,2Xix 2 + AH, 
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Fig. 1 




Fig. 2 



Fig. 1 - "Star-triangle" relation to reduce a three-connected spin xo- 
obtained in Eq. 



The new bonds below are 



Fig. 2 - Log-log plot of the variance cr(AE) of the defect energy as a function of systems size L for 
various bond fractions p > p c in d — 3 (left) and 6 (right). At small p, a(AE) drops to zero rapidly 
for increasing L, but turns around and rises for larger p, indicative of a nontrivial ordered state at 
low T. The longest transients in a(AE) mark p* , suggesting p^=3 = 0.272(1) and p^=6 = 0.0950(3). 



Jl,2 = § + Ji,2\ - \Ji,l - Ji,2\) > AH =^ (I J M + J iM + \Ji,l - Ji,2\) ■ ( 3 ) 

Rule V: A three-connected spin i can be reduced via a "star-triangle" relation, see Fig. ^ 



J M Xi X x + J i>2 X t X 2 + Ji,3 Xi x 3 < J 1)2 Xi x 2 + Ji, 3 II ^3 + ^2,3 ^2 ^3 + AiT, 
Ji, 2 = -A- B + C + D, J 1)3 =A-B + C-A J 2 , 3 = -A + B + C-Z? 
Ai? = A + B + C + A A=i|J a - J 4i2 + J ii3 |, 

^ = I 1^,1 — Ji,2 ~ Ji,3\ i C = 4 + ^i,2 + «^»,3 ] j B = j | Jj j + J, 



(4) 



8,2 



i,3 



The bounds in Eqs. 1314(1 become exact when the remaining graph takes on its ground state. 
Reducing higher-connected spins leads to (hypcr-)bonds between multiple spins, unlike Eq. 

After a recursive application of these rules, the original lattice graph is either completely 
reduced (which is almost always the case for p < p c ), in which case H D provides the exact 
ground state energy already, or we are left with a highly reduced, compact graph in which no 
spin has less than four connections. We obtain the ground state of the reduced graph with 
the extremal optimization heuristic [22], which together with H provides a very accurate 
approximation to the ground state energy of the original diluted lattice instance. 

Numerical Results. - In Ref. [20] it was shown that spin glasses with a discrete bond 
distribution on diluted lattices may possess a distinct critical point p* in their bond density, 
which is related to the (purely topological) percolation threshold p c of the lattice and the 
distribution of the bond weights P(J). Clearly, no long-range correlated state can arise below 
p c . A critical point distinct from percolation, p* > p c , emerges when such a correlated state 
above p c remains suppressed due to collaborative effects between bonds [20] (see Rule III). 
Thus, to observe any glassy properties on a dilute lattice, we have to ascertain p > p*. In 
Ref. [18], we were able to locate p* for the Migdal-Kadanoff lattice, by using the defect energy 
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Fig. 3 - Log-log plot of the width a of the defect energy distribution as a function of system size L 
for d — 3 (top) and 6 (bottom). The data is grouped into sets (connected by lines) parameterized 
by the bond density p. Most sets show a distinct scaling regime as in Eq. Q for a range of L above 
finite scaling corrections but below failing heuristic accuracy. 



scaling from Eq. QJ: For all p > p* the stiffness exponent y eventually took on its p = 1 value, 
while for any p < p* defect energies diminished rapidly for increasing L. 

In each dimension (see Ref. [9]), we have run the above algorithm on a large number of 
graphs (about 10 5 — 10 6 for each L and p) for p increasing from p c in small steps. For each 
given p, L increased until it was clear that a(AE) would either drop or rise for good. In this 
way, we bracket-in p*, as shown in Figs. |21 for d — 3 and 6. Tab. [I] summarizes those results. 

While the value of p* is distribution-dependent, we merely had to establish a bond fraction 
beyond which we would expect Eq. Q to hold. Then we can conduct numerical experiments 
to extract the asymptotic scaling of a(AE) for conveniently chosen p. For each choice of L 
and p, we have sampled the defect energy distribution with at least N > 10 5 instances, then 
determined its variance a(AE). Throughout, the distribution of AE is approximately normal, 
giving rise to an error bar oc 1/yN. In Fig. El we plot all the data for each dimension simply 
according to Eq. Q), on a logarithmic scale. For most sets of graphs, a scaling regime (linear 
on this scale) is visible. Yet, various deviations from scaling can be observed. Clearly, each 
sequence of points should exhibit some form of finite size corrections to scaling for smaller 
L. For large L, the inability to determine defect energies correctly, will inevitably lead to 
a systematic bias in a. Some data sets did not exhibit any discernible scaling regime, most 
notably our set for the undiluted lattice in d = 3. This resembles the observation of Refs. [23], 
which found long transients in similar studies on undiluted d — 2 or Migdal-Kadanoff lattices. 

To obtain an optimal scaling collapse of the data, we focus on the data inside the scaling 
regime for each set. To this end, we chose for each data set a lower cut in L by inspection. An 
appropriate high-end cut is introduced by eliminating all data points for which the remainder 
graph had a typical size of > 700 spins; at that point the EO heuristic (within the supplied 
runtime) seems to fail in determining defect energies with sufficient accuracy. All the remaining 
data points for L and p are fitted to a four-parameter scaling form, 

a(AE)^y ^L( P -p*) v 'Y (L-foo), (5) 

suggested by Ref. [20]. Unfortunately, we have a-priori no knowledge of corrections for small 
systems, making the low-L cut on the data a necessity. The fitted values for the fitting 
constants 3>o, v * ' , and y are listed in Tab.[l] Using the parameters of this fit, we re-plot only 
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L(p-p*) v L(p-p) v ' 



Fig. 4 - Scaling plot of the data from Figs|!|]for a, fitted to the functional form in Eq. JSJ as a function 
of the scaling variable x = L(p — p*)" for d — 3 (top) and 6 (bottom). Data that fell outside of 
the scaling regime in each set of Figs. [3] was cut. The straight line in each case represents the fit 
according to Eq. @ for the optimal collapse of the data, which provides an accurate determination 
of the stiffness exponent y in each dimension. For d — 3, we have also included the data for p — 1, 
which does not appear to have a scaling regime. 

the data from the scaling regime in Figs.0J 

We have also considered the case of d = 7. The implementation of our algorithm is not yet 
capable of handling lattice graphs with much more than 10 6 spins, limiting attainable sizes 
to L < 8, which was insufficient to determine p* directly (as in Fig- as wei l as the low end 
of the scaling variable x — L{p — p* ) u . Yet, we have computed about as many data points as 
for d = 6 (see "DoF" in Tab.HJ), as shown in the scaling plot in Fig. [5J Solid scaling extends 
over nearly half a decade, yielding j/7 = 1.24(5). Assuming that the data presented in Fig. [S] 
is asymptotic, any value for 2/7 outside of its error bars has a "goodness of fit" Q [24] below 
0.001, and Q < 10~ 100 for 1 or 3/2. Generally, yd appears to be rising with d above the upper 
critical dimension, but with values and at a rate below replica theory predictions [10]. 

* * * 
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Table I - List of the values for the critical bond- density p* (from Fig. \^for d = 3,6), and the fitted 
values according to Eq. of the correlation-length exponents v* , the surface tensions 3^o, and the 
stiffness exponents y, followed by the "goodness of fit," Q [24], and the degrees of freedom (DoF) used 
in the fit. Also listed are values for d — 4 and 5 from Ref. [9] and the fitted values for d = 7. 



d 


P* 


V* 


y 


y 


Q(DoF) 


3 


0.272(1) 


1.17 


2.37 


0.239 


1.00(92) 


4 


0.1655(5) 


0.60 


2.43 


0.610 


0.00(47) 


5 


0.1206(2) 


0.50 


3.05 


0.876 


0.86(48) 


6 


0.0950(3) 


0.44 


3.87 


1.103 


0.02(46) 


7 


0.080(1) 


0.40 


5.18 


1.244 


1.00(55) 
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Fig. 5 - Scaling plot of the data for a in d — 7 dimensions, fitted to the functional form in Eq. JjjJ 
as a function of the scaling variable x = L(p — p*) v . As before, data outside of the scaling regime, 
although plotted here, was disregarded in the fit. The straight line represents the fit according to 
Eq. @ for the optimal collapse of the data. It seems to exclude a simple (half-)integer value for 1/7, 
such as unity or 3/2 (dashed lines). 
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